MAE augmentation (primary and secondary outcome definitions)

library(tidyverse)
library(magrittr)
library(gt)
library(patchwork)
library(SummarizedExperiment)
library(tidySummarizedExperiment)
library(MultiAssayExperiment)

tmp <- fs::dir_map("R scripts/", source)
tmp <- fs::dir_map("../VIBRANT-99-utils/R/", source)
rm(tmp)
theme_set(theme_light())

Loading the MultiAssayExperiment object

mae_files <- 
  fs::dir_ls(
    get_02_output_dir(),
    regexp = "mae_full_.*\\.rds$"
  ) |> 
  sort(decreasing = TRUE) |>
  magrittr::extract(1)

mae <- readRDS(mae_files)

rm(mae_files)

Primary outcomes

Primary outcome definition

The primary outcome is “colonization with any of the L. crispatus strains contained in the LBP by 5 weeks of follow-up as assessed by metagenomic sequencing of the vaginal microbial community with detection of any one of LBP strains at >5% relative abundance or a combination of the strains accounting for >10% relative abundance by metagenomics.”

By metagenomics

“LBP colonization” at each visit

We compute, at each visit and for each participant, the total relative abundance of all LBP strains or the maximum relative abundance of any LBP strain.

If the total relative abundance of all LBP strains is larger than 0.1 or the maximum relative abundance of any LBP strain is larger than 0.05, we consider that the LBP achieved colonization_mg at that visit.

colonization_LBP_mg <- 
  mae[["mg"]] |> 
  as_tibble() |>
  left_join(mae@colData |> as_tibble()) |> 
  dplyr::filter(sample_type == "Clinical sample", !poor_quality_mg_data) |> 
  dplyr::filter(!is.na(LBP)) |> 
  group_by(.sample) |> 
  summarize(
    n = n(),
    total_LBP_rel_ab_at = sum(rel_ab),
    max_LBP_rel_ab_at = max(rel_ab),
    n_detected_strains_at = sum(rel_ab > 0),
    n_detected_LC106_strains_at = sum(rel_ab[LBP == "LC-106 & LC-115"] > 0),
    n_detected_LC115_only_strains_at = sum(rel_ab[LBP == "LC-115"] > 0)
  ) |> 
  mutate(
    colonized_LBP_at_mg = (total_LBP_rel_ab_at > 0.1) | (max_LBP_rel_ab_at > 0.05)
  ) 

colonization_LBP_mg <-  
  colonization_LBP_mg |> 
  dplyr::left_join(
    mae@colData |> as_tibble() |> 
      select(uid, pid, visit_code, randomized, arm, location) |> 
      dplyr::rename(.sample = uid), 
    by = join_by(.sample)
  ) |> 
  arrange(pid, visit_code) 

colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(    
    strip.text.y = element_text(angle = 0)
  )

“LBP colonization” by each visit

From the colonization_mg status at each visit, we can compute the colonization_mg status by each visit.

A participant is considered to have colonized by a visit if they have been colonized at that visit or any previous visit, starting from their post-product visit (“1200” for all groups).

Since we have some missing visits, we impute these as “not colonized”

colonization_LBP_mg <- 
  colonization_LBP_mg |> 
  select(-c(arm, randomized, location)) |> 
  dplyr::full_join(
    colonization_LBP_mg |> select(pid, arm, randomized, location) |> distinct() |> 
      expand_grid(
        visit_code = unique(colonization_LBP_mg$visit_code) |> sort()
      ),
    by = join_by(pid, visit_code)
  )

colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

colonization_LBP_mg <- 
  colonization_LBP_mg |> 
  mutate(colonized_LBP_at_mg = colonized_LBP_at_mg |> replace_na(FALSE)) |>
  arrange(pid, visit_code) |>
  group_by(pid) |> 
  mutate(
    tmp = ifelse(as.numeric(visit_code) < 1200, FALSE, colonized_LBP_at_mg),
    colonized_LBP_by_mg = cummax(tmp) |> as.logical()
  ) |> 
  ungroup() |> 
  select(-tmp)

colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_by_mg) +
  geom_tile() +
  annotate(
    geom = "rect", 
    xmin = 10.5, xmax = 11.5, ymin = -Inf, ymax = Inf,
    col = "black", alpha = 0, linewidth = 1
  ) +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

“L. crispatus colonization” at each visit

We compute, at each visit and for each participant, the total relative abundance of L. crispatus.

If the total relative abundance of L. crispatus strains is larger than 50%, we consider that the participant’s VM was colonized by L. crispatus at that visit.

colonization_crispatus_mg <- 
  mae[["mg"]] |> 
  as_tibble() |>
  left_join(mae@colData |> as_tibble()) |>
  dplyr::filter(sample_type == "Clinical sample", !poor_quality_mg_data) |> 
  dplyr::filter(species == "crispatus") |> 
  group_by(.sample) |> 
  summarize(
    total_crispatus_rel_ab_at = sum(rel_ab),
  ) |> 
  mutate(
    crispatus_dominance_at_mg = (total_crispatus_rel_ab_at >= 0.5)
  ) 

colonization_crispatus_mg <-  
  colonization_crispatus_mg |> 
  dplyr::left_join(
    mae@colData |> as_tibble() |> 
      select(uid, pid, visit_code, randomized, arm, location) |> 
      dplyr::rename(.sample = uid), 
    by = join_by(.sample)
  ) |> 
  arrange(pid, visit_code) 

colonization_crispatus_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = crispatus_dominance_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(    
    strip.text.y = element_text(angle = 0)
  )

“Crispatus colonization” by each visit

From the colonization_mg status at each visit, we can compute the colonization_mg status by each visit.

A participant is considered to have colonized by a visit if they have been colonized at that visit or any previous visit, starting from their post-product visit (“1200” for all groups).

Since we have some missing visits, we impute these as “not colonized”

colonization_crispatus_mg <- 
  colonization_crispatus_mg |> 
  select(-c(arm, randomized, location)) |> 
  dplyr::full_join(
    colonization_crispatus_mg |> select(pid, arm, randomized, location) |> distinct() |> 
      expand_grid(
        visit_code = unique(colonization_crispatus_mg$visit_code) |> sort()
      ),
    by = join_by(pid, visit_code)
  )

colonization_crispatus_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = crispatus_dominance_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

colonization_crispatus_mg <- 
  colonization_crispatus_mg |> 
  arrange(pid, visit_code) |>
  mutate(crispatus_dominance_at_mg = crispatus_dominance_at_mg |> replace_na(FALSE)) |>
  group_by(pid) |> 
  mutate(
    tmp = ifelse(as.numeric(visit_code) < 1200, FALSE, crispatus_dominance_at_mg),
    crispatus_dominance_by_mg = cummax(tmp) |> as.logical()
  ) |> 
  ungroup() #|> 
  #select(-tmp)

colonization_crispatus_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = crispatus_dominance_by_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

LBP colonization VS crispatus dominance

p1 <- colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

p2 <- colonization_crispatus_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = crispatus_dominance_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

 p1 + p2 + 
    plot_annotation(
      title = str_c("LBP colonization VS crispatus dominance")
    ) +
    plot_layout(ncol = 2) &
    theme(
      legend.justification = "left"
    )

p1 <- colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_by_mg) +
  geom_tile() +
  annotate(
    geom = "rect", 
    xmin = 10.5, xmax = 11.5, ymin = -Inf, ymax = Inf,
    col = "black", alpha = 0, linewidth = 1
  ) +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

p2 <- colonization_crispatus_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = crispatus_dominance_by_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

p1 + p2 + 
    plot_annotation(
      title = str_c("LBP colonization VS crispatus dominance")
    ) +
    plot_layout(ncol = 2) &
    theme(
      legend.justification = "left"
    )

By qPCR

“LBP colonization” at each visit

We compute, at each visit and for each participant, the relative abundance of each LBP strain as estimated by qPCR by dividing the copies per swab for that strain by the copies per swab for the 16S rRNA gene.

Then, just as we did for the metagenomics, we compute the total relative abundance of all LBP strains or the maximum relative abundance of any LBP strain.

If the total relative abundance of all LBP strains is larger than 0.1 or the maximum relative abundance of any LBP strain is larger than 0.05, we consider that the LBP achieved colonization at that visit by qPCR.

We also consider an additional secondary outcome, which is that if at least 2 strains are detected at >0, we consider that the participant is colonized by LBP at that visit by qPCR.

colonization_qPCR <- 
  mae[["qPCR"]] |> 
  as_tibble() |>
  left_join(mae@colData |> as_tibble()) |>
  dplyr::filter(sample_type == "Clinical sample") |> 
  dplyr::filter(!is.na(LBP)) |> 
  # we add the 16S copies per swab to compute the relative abundance
  left_join(
    mae[["qPCR"]] |> 
      as_tibble() |>
      left_join(mae@colData |> as_tibble()) |>
      dplyr::filter(sample_type == "Clinical sample") |> 
      dplyr::filter(is.na(LBP), .feature == "16S") |> 
      select(.sample, copies_per_swab_med) |>
      dplyr::rename(copies_per_swab_16S = copies_per_swab_med),
    by = join_by(.sample)
  ) |> 
  mutate(rel_ab = copies_per_swab_med / copies_per_swab_16S) |> 
  select(.feature, .sample, rel_ab, copies_per_swab_med, copies_per_swab_16S, everything()) |> 
  group_by(.sample) |> 
  summarize(
    total_LBP_rel_ab_at = sum(rel_ab),
    max_LBP_rel_ab_at = max(rel_ab),
    n_detected_strains_at = sum(copies_per_swab_med > 0),
    n_detected_LC106_strains_at = sum(copies_per_swab_med[LBP == "LC-106 & LC-115"] > 0),
    n_detected_LC115_only_strains_at = sum(copies_per_swab_med[LBP == "LC-115"] > 0)
  ) |> 
  mutate(
    colonized_LBP_at_qpcr = (total_LBP_rel_ab_at > 0.1) | (max_LBP_rel_ab_at > 0.05),
    LBP_detected_at_qpcr = (n_detected_strains_at > 2)
  ) 

colonization_qPCR <-  
  colonization_qPCR |> 
  dplyr::left_join(
    mae@colData |> as_tibble() |> 
      select(uid, pid, visit_code, randomized, arm, location) |> 
      dplyr::rename(.sample = uid), 
    by = join_by(.sample)
  ) |> 
  arrange(pid, visit_code) 

colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_qpcr) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(    
    strip.text.y = element_text(angle = 0)
  )

colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = LBP_detected_at_qpcr) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(    
    strip.text.y = element_text(angle = 0)
  )

“LBP colonization” by each visit

From the colonization status at each visit, we can compute the colonization status by each visit.

A participant is considered to have colonized by a visit if they have been colonized at that visit or any previous visit, starting from their post-product visit (“1200” for all groups).

Since we have some missing visits, we impute these as “not colonized”

colonization_qPCR <- 
  colonization_qPCR |> 
  select(-c(arm, randomized, location)) |> 
  dplyr::full_join(
    colonization_qPCR |> select(pid, arm, randomized, location) |> distinct() |> 
      expand_grid(
        visit_code = unique(colonization_qPCR$visit_code) |> sort()
      ),
    by = join_by(pid, visit_code)
  )

colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_qpcr) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(    
    strip.text.y = element_text(angle = 0)
  )

colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = LBP_detected_at_qpcr) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

colonization_qPCR <- 
  colonization_qPCR |> 
  arrange(pid, visit_code) |>
  mutate(
    colonized_LBP_at_qpcr = colonized_LBP_at_qpcr |> replace_na(FALSE),
    LBP_detected_at_qpcr = LBP_detected_at_qpcr |> replace_na(FALSE)
    ) |>
  group_by(pid) |> 
  mutate(
    tmp = ifelse(as.numeric(visit_code) < 1200, FALSE, colonized_LBP_at_qpcr),
    colonized_LBP_by_qpcr = cummax(tmp) |> as.logical(),
    tmp = ifelse(as.numeric(visit_code) < 1200, FALSE, LBP_detected_at_qpcr),
    LBP_detected_by_qpcr = cummax(tmp) |> as.logical()
  ) |> 
  ungroup() |> 
  select(-tmp)

colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_by_qpcr) +
  geom_tile() +
  annotate(
    geom = "rect", 
    xmin = 10.5, xmax = 11.5, ymin = -Inf, ymax = Inf,
    col = "black", alpha = 0, linewidth = 1
  ) +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = LBP_detected_by_qpcr) +
  geom_tile() +
  annotate(
    geom = "rect", 
    xmin = 10.5, xmax = 11.5, ymin = -Inf, ymax = Inf,
    col = "black", alpha = 0, linewidth = 1
  ) +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

Colonization qPCR VS metagenomics

p1 <- colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_mg) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

p2 <- colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_at_qpcr) +
  geom_tile() +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
    ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

p1 + p2 + 
    plot_annotation(
      title = str_c("LBP colonization *at* metagenomics VS qPCR")
    ) +
    plot_layout(ncol = 2) &
    theme(
      legend.justification = "left"
    )

p1 <- colonization_LBP_mg |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_by_mg) +
  geom_tile() +
  annotate(
    geom = "rect", 
    xmin = 10.5, xmax = 11.5, ymin = -Inf, ymax = Inf,
    col = "black", alpha = 0, linewidth = 1
  ) +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )


p2 <- colonization_qPCR |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = colonized_LBP_by_qpcr) +
  geom_tile() +
  annotate(
    geom = "rect", 
    xmin = 10.5, xmax = 11.5, ymin = -Inf, ymax = Inf,
    col = "black", alpha = 0, linewidth = 1
  ) +
  facet_grid(
    randomized + arm ~ ., 
    scales = "free_y", space = "free_y", labeller = label_both
  ) +
  theme(
    strip.text.y = element_text(angle = 0)
  )

p1 + p2 + 
    plot_annotation(
      title = str_c("LBP colonization by metagenomics VS qPCR")
    ) +
    plot_layout(ncol = 2) &
    theme(
      legend.justification = "left"
    )

Adding primary_outcomes assay to MAE

se_primary_col_LBP_mg <-
  SummarizedExperiment(
    assays = list(
      outcome = 
        colonization_LBP_mg |> 
        # filter(!is.na(pid)) |> 
        # mutate(.sample = str_c(pid, "_", visit_code)) |> 
        dplyr::filter(!is.na(.sample)) |> 
        select(.sample, colonized_LBP_at_mg, colonized_LBP_by_mg) |> 
        as.data.frame() |> 
        column_to_rownames(".sample") |> t()
  ) 
)

se_primary_col_crisp_mg <-
  SummarizedExperiment(
    assays = list(
      outcome = 
        colonization_crispatus_mg |> 
        dplyr::filter(!is.na(.sample)) |> 
        select(.sample, crispatus_dominance_at_mg, crispatus_dominance_by_mg) |> 
        as.data.frame() |> 
        column_to_rownames(".sample") |> t()
  ) 
)

se_primary_col_LBP_qPCR <-
  SummarizedExperiment(
    assays = list(
      outcome = 
        colonization_qPCR |> 
        dplyr::filter(!is.na(.sample)) |> 
        select(.sample, colonized_LBP_at_qpcr, colonized_LBP_by_qpcr, LBP_detected_at_qpcr, LBP_detected_by_qpcr) |> 
        as.data.frame() |> 
        column_to_rownames(".sample") |> t()
    )
  ) 


mae <- c(
  mae, 
  list(col_LBP_mg = se_primary_col_LBP_mg, 
       col_crispatus_mg = se_primary_col_crisp_mg,
       col_LBP_qPCR = se_primary_col_LBP_qPCR)
)

Saving the MultiAssayExperiment objects

We save the “master” MAE (with all assays); but if needed for publication, we can subset the MAE to only include the assays/data of interest.

saveRDS(
  mae, 
  str_c(
    get_03_output_dir(),  
    "mae_full_", today() |> str_remove_all("-"), ".rds"
  )
)